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Anharmonic properties from a generalized third order ab initio approach: theory and 

apphcations to graphite and graphene 

Lorenzo PaulattoQ Francesco Mauri, and Michele Lazzeri 
IMPMC, Universite Pierre et Marie Curie, CNRS, 4 place Jussieu, F-75005 Paris, France 

(Dated: April 10, 2013) 

We have implemented a generic method, based on the 2n+l theorem within density functional 
perturbation theory, to calculate the anharmonic scattering coefficients among three phonons with 
arbitrary wavevectors. The method is used to study the phonon broadening in graphite and graphene 
mono- and bi-layer. The broadening of the high-energy optical branches is highly nonuniform and 
presents a series of sudden steps and spikes. At finite temperature, the two linearly dispersive 
acoustic branches TA and LA of graphene have nonzero broadening for small wavevectors. The 
broadening in graphite and bi-layer graphene is, overall, very similar to the graphene one, the most 
remarkable feature being the broadening of the quasi acoustical ZO' branch. Finally, we study 
the intrinsic anharmonic contribution to the thermal conductivity of the three systems, within the 
single mode relaxation time approximation. We find the conductance to be in good agreement with 
experimental data for the out-of-plane direction but to underestimate it by a factor 2 in-plane. 

PACS numbers: 63.20.kg,63.20.dk,63.22.Rc,65.80.Ck 



I. INTRODUCTION 

Thermal transport is currently attracting much atten- 
tion; the main applications of interest are materials for 
thermoelectric energy conversion^ and materials used for 
thermal dissipation in microelectronics^. While, in the 
first case, the goal is to engineer the smallest possible 
thermal conduction, in the second case a good thermal 
conduction is required. In general, our understanding of 
thermal properties is heavily based on theoretical model- 
ing and the use of precise and reliable approaches, such as 
the ab initio computational methods, is highly desirable. 

The presence of a temperature gradient in a solid in- 
duces a heat flux. The heat carriers can be lattice vibra- 
tions (phonons) or electronic excitations. In general, lat- 
tice conduction is the dominant mechanism in the pres- 
ence of an electronic gap (semiconductors and insulators) 
or when the gap is zero but the density of electronic states 
at the Fermi level is small (semi-metals) . Lattice thermal 
resistance is then dictated by phonon scattering, which 
can be induced by extrinsic mechanisms (isotopic disor- 
der, structural defects, finite-size of the crystals, etc.) 
or by intrinsic ones (anharmonic phonon-phonon scat- 
tering) . Determining the intrinsic anharmonic scattering 
is in itself a very complex task; it has been attempted 
ab initio, within density fun ctio nal theory (DFT), us- 
ing finite difference derivatiorP'^ or molecular dynamics 
techniquea^ or from linear response theory.EH^l 

In a crystal, the intrinsic lattice thermal conduction 
can be obtained by knowing harmonic phonon energies 
and anharmonic phonon-phonon scattering coefficients. 
Harmonic phonon energies are determined by the sec- 
ond order derivative of the system total energy, with 
respect to atomic displacements. This second deriva- 
tive can be efficiently calculated ab initio by using den- 
sity functional perturbation theory (DFPT)p21^ which al- 
lows the determination of the phonon dynamical ma- 
trix for an arbitrary q wavevector in the Brillouin zone. 



DFPT is implemented in the Quantum ESPRESSO 
package^^, within the plane-waves and pseudopotential 
approaches. The anharmonic scattering coefficients can 
be determined by the third order derivative of the en- 
ergy with respect to three phonon perturbations, char- 
acterized by the wavevectors q, q', q". For the thermal 
transport problem, it is necessary to know these deriva- 
tives with respect to three arbitrary wavevectors (pos- 
sibly q 7^ 0, q' ^ 0, q" ^ 0), with the only condition 
q-l-q' + q" ~ G, where G is a reciprocal lattice vector. In 
principle, these coefficients can be obtained within DFPT 
by using the, so called, "2n-|-l" theorem as formulated by 
Ref. [121 This theorem allows us to access the 3rd deriva- 
tive of the total energy by using only the 1st derivative of 
ground state density and wavefunctions; contrary to the 
finite differences approach, we do not have to perform 
expensive supercell calculations. 

The first implementations of the "2n-|-l" approach 
were limited to the scattering of one zero-momentum 
phonon towards two phonons with opposite arbitrary mo- 
menta (0, -q, q). Ref . U^ implemented this approach for 
insulating and semiconducting materials. Later, Ref. [TJI 
generalized the approach to metals and zero-gap materi- 
als. The (0, -q, q) anharmonic coefficients can be com- 
puted within Quantum ESPRESSO ^^, using the d3 
code which was developed in ReL 14, These coefficients 
can be easily used to compute the anharmonic broaden- 
ing of a q = phonon (see e.g. Ref. 15). By using a 
super-cell approach one can also compute the broaden- 
ing of phonons having q commensurate with the super 
cell. However, the super-cell approach (which was used 
in Ref. 'W, 'W, and T5|) can be computationally very de- 
manding. Recently, Ref. [12] further extend the method 
to three arbitrary phonons, (q, q', q"), although only for 
insulators/semiconductors. This was done within a pro- 
prietary non-publicly-available software. 

We have developed and implemented an extension of 
the d3 code of the Quantum ESPRESSO software 



to compute, within the DFPT "2n+l" approach, the 
three-phonons anharmonic coefficients for three arbitrary 
wavevectors (q, q', q") for insulators/semiconductors and 
also for metaUic or zero gap systems. The coefficients 
thus obtained can be used in a straightforward way to 
compute the anharmonic broadening of a phonon with 
an arbitrary wavevector q and the intrinsic thermal con- 
ductivity within the single-mode relaxation time approx- 
imation (SMA). fiZl xhe ffist applications of the method 
are devoted to graphite, graphene and graphene bi-layer. 
Indeed, the thermal properti es of these systems have at- 
tracted significant attention,E^Iini being sp^ carbon sys- 
tems excellent thermal conductors. 

Sec. [n] describes the method. Sec. |III| reports the re- 
sults and the discussion. Conclusions are summarized in 
SecHVl 



II. METHOD 



where 



In Sec. II A we provide the expressions for the phonon 
anharmonic broadening and for the phonon thermal con- 
duction. In Sec. |IIB| the method is briefly described and 
we provide the relevant computational details. A more 
detailed description of the method is reported in the Ap- 
pendix. 



A. Anharmonic decay and thermal transport 

Let us consider the total energy for a crystal 
f*°*({iiR s Q,}), where wr,s,q is the displacement from the 
equilibrium position of the s atom along the a Cartesian 
coordinate in the unit cell identified by the lattice vector 
R. We define 
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where the sum is performed on the lattice vectors {R} 
and N is the number of cells involved in the summation. 
We define the dynamical matrix 
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the angular frequency ojq j of a phonon with wavevector 
q and branch index j is obtained by solving 
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V^'^' has the dimension of an energy and does not depend 
on N, while Xqj is adimensional. Because of the transla- 
tional symmetry of the crystal, the coefficients V^^' from 
Eq. |4] are 7^ only when q -I- q' -I- q" = G, where G is 
any reciprocal lattice vector. 

With these definitions, the lifetime due to anharmonic 
phonon-phonon interaction, Tq,, and the correspond- 
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Where T is the temperature, Tiqj is the Bose-Einstein 
statistics occupation of phonon (qj'), and 5{x) is the 
Dirac distribution. The sum is performed over a suf- 
ficiently fine grid of Nq q-points in the Brillouin zone 
(BZ) and q" — — q — q'. Tqj and 7qj depend on T only 
through the phonon occupations n. 

The r.h.s. of Eq. [6]is usually interpreted as the sum of 
scattering processes in which a phonon of wavevector q 
decays into two phonons — q', — q", (third line of Eq. [6]) 
or in which the phonon q coalesces with — q' and emits 
— q" (fourth line of Eq. [6]). The energy conservation of 
the processes are guaranteed by the Dirac delta. One can 
also distinguish between Normal and Umklapp processes: 
by choosing q and — q' such that they belong to the first 
BZ, the scattering is Normal when also q" = — q — q' 
belongs to the first BZ; on the contrary, when q" does 
not belong to the first BZ, the scattering is Umklapp. 

By knowing the anharmonic scattering coefficients, 
Eq. |4] one can determine the lattice thermal conductivity 
within the framework of the Boltzmann transport equa- 
tion (BTE) for phonons^'^. In general, an exact solution 
of the BTE is a difficult task; a commonly used approxi- 
mation to the problem is the so-c alled single mode relax- 
ation time approximation (SMAJ^EH. Within the SMA, 
the lattice thermal conductivity tensor becomes: 



where z are the orthogonal phonon eigenmodes normal- 
ized in the unit cell and nis is the atom mass. We define 
the three-phonon scattering coefficients as 
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Here, ft is the volume of the unit cell, K^ is the 
Boltzmann constant and c" ■ is the phonon group ve- 



locity of mode (qj) along Cartesian direction a: 



qj 



dujqj/{dqa). The SMA conductivity from Eq. 7] can be 
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obtained in a straightforward way once the an 
lifetimes t^j have been computed from Eq. 
is a 3 X 3 tensor which takes into account the possi- 
ble anisotropics and transversal conductance. However, 
in high-symmetry crystals, as graphene and graphite, 
the off-diagonal elements are zero, if two axes lie in 
the graphene plane. Moreover, in both graphene and 
graphite the two in-plane xx and yy components are iden- 
tical. The out-of-plane zz component is not defined in 
the bidimensional graphene systems, but it is meaningful 
in graphite. 

T he va lidity limits of the SMA are discussed in liter- 
ature^E3 Here, we wish to remind that, for a generic 
material, the SMA is expected to be valid (that is, to 
provide the correct solution to the BTE) at room condi- 
tions and to break down only at very small temperatures 
(see, e.g. Ref. O. 



B. Density-functional theory calculation 

Calculations of the phonon properties are done within 
density functional perturbation theory^" as implemented 
in Ref. JT., The third order coefficients defined in Eq. |4] 
are computed using a code which has been developed for 
the present work. This code has been written on the top 
of a previous less general implementation available within 
the Quantum ESPRESSO package: the d3 code, which 
was implemented in Ref. 13! The method is described in 
detail in Appendix |A] and [B| 

We use local-density approximation and the carbon 
atom is described by a norm-conserving pseudopotential 
which includes four electrons in valence. Plane waves 
kinetic energy cutoff is 90 Ry. For all the systems, the 
in-plane lattice parameter is a = 2.44 A, which is the the- 
oretical equilibrium value for graphite. For graphite, we 
use c/a = 2.664. This value, which is only slightly differ- 
ent from the experimental value c/a — 2.727, is chosen 
phenomenological to accurately reproduce the low fre- 
quency phonon dispersion along the F-A direction (see 
the discussion below). For graphene, the periodic repli- 
cas of the planes are spaced along the z direction with 
7 A of vacuum. The two layers of the graphene bilayer are 
spaced with the inter-planar distance of bulk graphite; 
periodic images are then spaced with 7 A of vacuum. 

The computational parameters are listed in Ap- 
pendix|D] We remind here that the electronic integration 
has to be done with a small value of smearing (and a con- 
sequent fine-grained k-point grid) due to the presence of 
a Kohn anomaly for the highest optical branch near K^'^ 
(usually called TO). The phonon frequencies Wqw and the 
third-order coefficients V^^\ used in Eqs.[6|and[7J are cal- 
culated in a slightly different way. On one hand, phonon 
energies are corrected using an ad hoc procedure (based 
on DFT-I-GW renormalization of the electron-phonon in- 
teraction as in Ref. 23 see Appendix [PJ). This correction 
affects only the TO branch, it does not touch the other 
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FIG. 1. (Color online) Graphite phonon dispersion. Sym- 
bols are measurements from Refs. [21] and 1221 Lines are cal- 
culations. In the upper panel, the TO optical branches are 
plotted two times. The solid (black) lines include GW cor- 
rections of the electron-phonon interaction. On the contrary 
the dashed (green) lines are done using standard DFT and 
they are shown only for comparison. Lower panel: Solid lines 
are done by using cja = 2.664; dashed lines are done with 
cja — 2.727 and they are shown only for comparison. In the 
inset we recall the high symmetry points naming convention. 
In both panels, the solid lines correspond to the calculations 
used throughout the paper. 



branches, and it provides better agreement with mea- 
surements. Fig. [l| On the other hand, the third-order 
coefficients are computed within standard (less precise) 
DFT. The using of these two different procedures for the 
LOqj and l^'^-* calculations is not consistent. However, 
this should not affect the results in a major way: phonon 
broadening results from a sum over different processes 
which are selected by energy conservation enforced by 
the two Dirac 6 in Eq. [6j The intensity of the processes 
is then proportional to the square of the V^"^^ coefficients. 



Consequently, the computational accuracy of ujqj and 
that of the T^'^^ coefficients affect the result in a very 
different way. An error in the phonon dispersion can af- 
fect the lifetime in a not predictable way and, thus, a 
special care should be taken into finding the best possi- 
ble description of the phonon dispersion. The same care 
is not strictly necessary for the third order calculations. 
Fig. IT] compares measured with calculated phonon dis- 
persions for graphite. Notice that plain DFT calculations 



do not provide a satisfactory description of the highest 
optical TO branch near K, while DFT+GW ones do 
much better. The lower panel of Fig. [T] shows in detail 
the low frequency dispersion. This region is character- 
ized by the splitting of the acoustic phonon branches of 
the two graphene planes in the graphite unit cell. These 
branches are particularly sensitive to the actual value of 
cj a. In particular, in that region, by changing the lattice 
parameters from cja = 2.664 (which is the value used 
throughout the paper) to cja = 2.727 the value of the 
phonon branches change by almost 14%. 

Actual DFT calculations are done on a relatively 
coarse grid of q wavevectors, described in Appendix [D] 
The dynamical matrices and the third order coefficients, 
necessary to compute the broadening and the thermal 
conductivity (Eqs.p^and ^, are then obtained on a finer 
grid via the Fourier interpolation technique described in 
Appendix [cj Eqs. [6]and [7] are evaluated by performing 
the sum over a discrete grid of q points and by substi- 
tuting the 6{x) with a Gaussian function characterized 
by an artificial smearing x- This approximation is valid 
as long as x is smaller than the thermodynamic fluctua- 
tion, which is of order K^T . The grids and the x values 
are specified in Appendix [Pj He re we just remark that 
the results shown in Sects. iTlIAl and IIIIBI are obtained 
using a particularly fine-grained sampling. This is only 
necessary to produce the very sharp features which are 
present in the broadening of the higher optical bands or 
to produce the correct behavior of the broadening of the 
acoustic branches in the vicinity of T. Indeed, a much 
coarser grid is sufficient for most applications, included 
those described in Sec. IIII CI 



III. RESULTS AND DISCUSSION 

This section reports and describes the results. Sec- 
tion [III A| analyzes the anharmonic phonon broadening 
in the graphene monolayer. Special relevance is given to 
the three acoustic branches which are the most impor- 
tant for the thermal transport. Sec. |IIIB| analyzes the 
broadening in graphene bilayer and graphite. Sec. |III C| 
is dedicated to thermal transport. 



A. Graphene phonon broadening 

Fig. [2] shows the calculated graphene phonon disper- 
sion, the respective anharmonic broadening and the vi- 
brational density of states (VDOS). The branches are 
labeled in the usual way.l^ There are three acoustic 
branches (ZA, TA, LA) and three optical branches (ZO, 
TO, LO). ZA and ZO correspond to an atomic motion 
perpendicular to the graphene plane {z direction), all the 
other branches are polarized parallel to the plane. In the 
vicinity of F, TA and TO are quasi transverse, while LA 
and LO are quasi longitudinal. In the following, these 
labels will be used to classify the branches all along the 



high symmetry lines (as in Fig. [2]), although this dis- 
tinction is not meaningful for an arbitrary wavevector in 
the Brillouin zone (BZ). Because of symmetry, the modes 
perpendicular polarized (ZA and ZO) are separated from 
the others all over the BZ. Moreover, the TA branch is 
always well separated from the other parallel polarized 
branches (labeled as Fh)- In Fig. [2] we can, thus, sepa- 
rate the VDOS in three distinct components labeled as 
Z, TA, and Pr- The two dimensional character of the 
phonon dispersion is associated with some specific fea- 
tures. The ZA branch is quadratic near T and, thus, in 
the limit w — >■ the VDOS does not go to zero (Fig. |2|. 
The presence of a local maximum in the phonon disper- 
sion (as the one at 1008 cm~^ for the TA branch near K 
or the one at 904 cm~^ for the ZO one near F) is associ- 
ated with a step in the VDOS. The presence of a saddle 
point in the dispersion (as those at 477 cm^^, 631 cm~^, 
643 cm^^, and 1432 cm^^ at the M point) is associated 
with a sharp peak in the VDOS. 

Fig.lSlreports in more detail the calculated anharmonic 
phonon broadening, along high symmetry lines, and its 
decomposition into the different allowed decay channels. 
For symmetry reasons, the z-polarized branche s can only 
decay toward one Z and one non-Z phonons.^SEI! xhe 
other bands can only decay towards two phonons which 
are either both or neither z-polarized. The two most 
striking features in Figs. [2]|3]are the small q behavior of 
the acoustic branches and the highly non uniform behav- 
ior of the broadening. 

First, we remind that in a three dimensional isotropic 
crystal, all the three acoustic branches are linearly dis- 
persive and one expects to observe for 5 —7> a vanishing 
broadening. On the contrary, at finite temperature, both 
TA and LA branches of the two dimensional graphene 
have a nonzero broadening in the g — > limit. This be- 
havior is due to a decay process in which a TA (or LA) 
phonon decays into two phonons both in the ZA branch. 
Fig. [3] This decay is entirely due to Normal scatter- 
ing. This can be seen in Fig. l4l where the broadening of 
the acoustic branches is decomposed into the two com- 
ponents which are due, respectively to Normal and Umk- 
lapp processes, as defined in Sec. |II A| Actually, one can 
demonstrate^ that, in general, when a linearly dispersive 
phonon decays into two quadratically dispersive phonons, 
the broadening is non vanishing in the g — )• limit be- 
cause of energy and momentum conservation. Moreover, 
the quadratically dispersive ZA branch has a broadening 
which is itself quadratic in q around F. The ZA broad- 
ening is due to a Normal decay in which the ZA phonon 
decays into one ZA phonon and one linearly dispersive, 
TA or LA, phonon. Again, one can demonstrate that, 
in general, when a quadratically dispersive phonon de- 
cays into one quadratically and one linearly dispersive 
phonons, the broadening vanishes quadratically in the 
q — > limit. We remark that, here, the anharmonic 
broadening has been computed by summing over an ex- 
tremely fine reciprocal-space grid (see AppendixO). This 
is necessary in order to reproduce correctly the anoma- 
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FIG. 2. (Color online) Calculated graphene phonon dispersion. Each phonon branch is represented with a variable-width filled 
band: the graphical width is equal to the respective anharmonic broadening at 300 K, expressed in cm~^ and magnified by a 
factor 100. The vibrational density of states (VDOS) is also shown, together with its decomposition over groups of disentangled 
branches labeled as Z (corresponding to the ZA and ZO branches), TA and Ph (corresponding to LA, TO, and LO). 




FIC. 3. (Color online) Graphene anharmonic phonon broad- 
ening (FWHM) at 300 K, for each phonon branch (labeled as 
in Fig. pi, along high symmetry lines. The total broadening 
(solid thick line) is decomposed depending on the character 
of the final states which are labeled as Z, TA, and Ph (see the 
text): e.g. TA-Ph corresponds to a decay involving one TA 
and one Ph phonon. 



lous behavior of the LA and TA broadening for small q. 
Far from F, the details of the broadening ean be correctly 
reproduced by using a much coarser grid and a larger x- 

The existence of a finite broadening at small q for the 
TA and LA acoustic branches is problematic. Indeed, 
the concept itself of phonon is meaningful only when 
w/7 > 1, being 7 the broadening (i.e. the inverse of 
the phonon lifetime). From the present calculations, the 
condition uj/j > 1 is satisfied for both the TA and LA 
branches for q > q, with q — 0.5 x 10~'*27r/ao, being qq 
the in-plane lattice spacing. Thus, for q < q, the TA or 
LA frequency can become smaller that the broadening. 
In this region the present treatment is, obviously, not 
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FIG. 4. (Color online) Graphene anharmonic phonon broad- 
ening for the three acoustic branches at 300 K (same data 
as in Fig. |3|: the broadening is decomposed in two compo- 
nents which are due, respectively, to Normal and Umklapp 
processes. 




FIG. 5. (Color online) Graphene anharmonic phonon broadening (FWHM, in cm ^) at 300 K, for the three acoustic branches, 
over all the Brillouin zone. 



valid (see the discussion in Ref. :9) and a proper treat- 
ment of the phenomenon is beyond the present scope. In 
practice, however, q <q represents a tiny portion of the 
Brillouin zone (the corresponding region in Figs. [2] p^ has 
width of the order of the thickness of the vertical line 
passing through T). As a consequence, we can assume 
that the properties obtained as a sum over the Brillouin 
zone (such as the thermal conductivity of Eq. n[ are not 
affected by a major error. 

Concerning the global appearance of Figs. [2j |3j the 
many sharp peaks in the broadening can be ascribed to 
different mechanisms. Those in the high energy part of 
the spectrum are, in general, associated with peaks in the 
VDOS: when one or both of the final states {i.e. of the 
states that meet the energy and momentum conservation 
requirements in Eq. pi) produce a peak in the VDOS, 
there the broadening will typically exhibit a peak. For 
example, the large scattering probability, predicted for 
the M point on the LO branch (1373 cm"^), corresponds 
to a decay toward a ZO phonon close to T (904 cm~^) 
and a ZA phonon close to M (477 cm^^). As the VDOS, 
Fig. [2j has maximum in both region, this transition is 
particularly favored. On the other hand, for q>0.66M 
(along the TM direction) or for q>0.62K (along TK), 
the LO broadening displays a sudden increase. This is 
because, for these wavevectors, the LO energy becomes 
small enough to activate the decay channel towards the 
ZO branch (At q-0.66M and q~0.62K the LO phonon 
decays into a ZO with the same wavevector and a ZO 
with q^ r). 

We remark that the presence of sharp peaks which are 
essentially determined by energy and momentum conser- 
vation in the decay process implies that even a small 
change in the phonon dispersion used in the calcula- 
tion could induce significant differences in the calculated 
broadening. 

Concerning the three acoustic branches the broaden- 
ing near T is almost entirely due to Normal scattering. 
Fig. W] The peaks which are observed at q>0.44 M and 
q>0.41 K for the LA branch, and at q>0.65 M and 
q>0.54 K for the TA one, are associated with the acti- 
vation of Umklapp scattering towards the ZA phonons. 
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FIG. 6. (Color online) Phonon mean free path for the three 
acoustic branches in graphene at 300 K. 



To have a more comprehensive view. Fig. [5] reports the 
broadening in the entire Brillouin zone. The TA and 
LA branches exhibit a feature-rich behavior in a wide re- 
gion, far from F, which roughly starts at about halfway 
to the first BZ edge. In this region, the anharmonic de- 
cay presents a component of Umklapp processes, which is 
absent in the vicinity of F, where the scattering is almost 
entirely Normal. On the other hand, the ZA broadening 
is relatively feature-less and isotropic; it is quadratic in q 
in the center of the first BZ then it saturates and remains 
roughly constant. 



1. Phonon mean free path 

An alternative way to represent the effect of the 
phonon broadening is to plot the single-phonon mean free 
path (MFP): 



-^qj — Tq 



where Cq, is the phonon group velocity and Tq 



(8) 

the 

phonon lifetime, from equation [6J In figure |6] we have 
plotted the MFP at 300 K for the three acoustic branches. 
The MFP for the TO and LO bands is of order 100 nm 
or smaller. We have verified that it does not get substan- 
tially higher at lower temperatures, except in the vicinity 
of F where it diverges at K. At room temperature, the 
MFP of the ZA bands is one order of magnitude larger. 
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FIG. 7. (Color online) Q and a parameters (defined in the text) for the graphene acoustic (left) and optical (right) 
branches. For a given phonon mode, the temperature dependence of the anharmonic broadening can be approximated by 
7(r) ~ aS coth(0 /T), where Q and a. are the parameters corresponding to that phonon. Q is a characteristic temperature 
and Q is the high temperature slope of y{T). 70 is the T = K broadening. 



i.e. of order 1 /im, in the center region of the Brillouin 
zone. Furthermore, it increases as 1/T when tempera- 
ture decreases. Obviously, especially at small temper- 
atures, when the intrinsic anharmonic MFP is too big, 
other effects (typically the scattering on the borders of 
the sample) become important and limit the value of the 
MFP. We remark that the MFP of acoustical phonons 
is only one order of magnitude smaller than the typi- 
cal dimensions of high-quality graphene samples. It is, 
also, definitely larger than the transverse dimension of 
graphene nano-ribbons. This result suggests that ballis- 
tic phonon-driven conductance could play a relevant role 
in this kind of systems. 



where 



'qj 



ttK 



B 



riv, 



E k 



' q'j'j" 

1 1 



(3) 

qi,q'j',q"j" 



■'q J ^q"j' 



6{UJ^-Ujq,j, -Wq"j") 



Oq,j> UJqnj,, 



Siujqj+UJq'j' -Wq"j") 



(10) 



does not depend on T. One can thus be tempted to ap- 
proximate the overall dependence on T of the broadening 

7 by 



2. Temperature dependence 

The intrinsic anharmonic broadening of a specific 
phonon (qj) has, in general, a typical dependence of the 
temperature T: It is almost constant below a certain 
characteristic temperature 9^, then it rapidly becomes 
linear in T. Such a behavior is reproduced by Eq. [6j A 
quadratic dependence on T can be observed experimen- 
tally only at relatively high T and it is due to terms of 
order higher than those included in Eq. [6J ESI 

From Eq. [6j one can check that 



lim 7q,(T)=aq,-r + 0(l/T), 



(9) 



7q, (r) =^ 7q,- (T) = aq, C Coth Mf 



(11) 



where coth is the hyperbolic cotangent function, 0^ = 
7qj(0)/Q;qj, and 7qj(0) is the T ~ broadening from 
Eq. |6] Indeed, 7 from Eq. [Tl] is almost constant for 
T <^ 9^ and tends to 7(0) for T ^ 0. Moreover, 
7(r) = aT -h 0(l/r) for T > 6lS. To check the va- 
lidity of this approximation (Eq. 11), we systematically 



computed the graphene broadening for different phonons 
in the temperature range between and 1500 K, using 
Eq. J6l These results are reasonably well reproduced by 
Eq. [11] with an error less than 5%. 

As a consequence, for a given phonon mode (qj), the 
knowledge of the two corresponding parameters 6 and a 
is enough to determine the overall temperature behavior 
of the broadening, by using Eq. 11 The two parameters 
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FIG. 8. (Color online) Graphite and graphene bilayer: calculated phonon dispersion widened by the anharmonic broadening 
(FWHM X 100) at 300 K. 



can be extracted from Fig. [7J for the graphene acoustic 
branches, along high symmetry lines. 

Finally, the broadening of the LA and TA branches 
can be fitted with an isotropic function of g = |q| and of 
T of the form: 



7(g,r) ==gBcoth 



qA 
T 



(12) 



By defining h^ — — ^ where ao is the cell parameter, 
we have: for the LA band, Sla — 4.58 cm~^/5o and 
Ai^A = 694 K/6o; for the TA band, B-^a = 0.805 cm^ V^o 
and A'YA = 241 K/6o. These fitted parameters reproduce 
the computed linewidth with an error of generally less 
than 10% for q < 0.40 bo (smaller for higher temperature 
and lower |q|). For the ZA band, 6^ — can be assumed, 
resulting in this simple fitting function: 



-f{q,T) = Bq^T . 



(13) 



Taking Bza = 25.9 x 10^'^ cm^^/^Q we reproduce 
the value of linewidth to 10% accuracy in the range 
q < 0.25 bo, except for systematically underestimat- 
ing it in the very small q region where it is negligible 
(7 < 10-5 cm-i). 



B. Graphite and bilayer graphene 

We now discuss the anharmonic broadening in graphite 
and graphene bilayer. Each of the six phonon branches 
of the graphene monolayer splits into two branches for 
both graphite and graphene bilayer. The three acous- 
tic branches of graphene (ZA, TA, LA) split into three 
acoustic (ZA, TA, LA) and three quasi-acoustic branches 
(ZO', TO', LO'). The quasi-acoustic branches are almost 
degenerate with their respective acoustic ones, except in 
the vicinity of the F-A line. The remaining six optical 
branches are pair- wise quasi-degenerate in the entire Bril- 
louin zone, and they will be referred in pairs simply as 
ZO, TO and LO or, in some cases as ZO^, ZO^, etc. This 
notation does not hold along the F— A line in graphite, as 
degeneracy changes. It is however still possible to name 
the branches by continuity below 600 cm"^. 



Fig. [8] shows a general view of the calculated phonon 
dispersion and broadening in bulk graphite and graphene 
bilayer. Fig. [9] and Fig. [TOJ compare the broadening of 
the acoustic and quasi-acoustic branches of, respectively, 
graphite and bilayer graphene with those of the single 
layer graphene. Fig. [Tl] shows in more detail the low fre- 
quency region. In the high energy part of the spectrum, 
graphite and graphene monolayer and bilayer are almost 
indistinguishable. Fig. IHl meaning that, also in graphite, 
the physics is ruled by the two dimensional character of 
the phonon dispersion. Some of the graphene sharp fea- 
tures are a slightly broader in graphite due to the out-of- 
plane phonon dispersion acting like an effective smearing. 
The most striking differences between three-dimensional 
bulk graphite and two-dimensional graphene monolayer 
and bilayer are associated with the acoustic and quasi- 
acoustic branches. 

First, we remind that, in graphene, the vibrational 
density of states (VDOS in Fig. [2]) has a finite constant 
value for energies approaching zero. This is because 
the ZA branch has a quadratic dispersion (not linear 
as usual) and the VDOS is calculated in a two dimen- 
sional Brillouin zone. On the contrary, in graphite, the 
VDOS goes to zero almost linearly for energies going to 
zero. Fig. [8j This happens in spite of the fact that in 
graphite, when the out-of-plane component q^ = 0, the 
ZA branch is not very different from the graphene one. In 
particular, on the scale of Fig. [8j the graphite ZA branch 
(for qz = 0) appears quadratic as the graphene one from 
Fig. [2] However, for graphite, the VDOS is calculated in 
a three dimensional BZ and the qz = phonons have an 
infinitesimal weight. Also, notice in the graphite VDOS, 
the presence of a peak at 132 cm~^ corresponding to the 
ZO' frequency at F. This peak and this phonon do not 
have a correspondence in graphene. Finally, from Fig. |8J 
the bilayer VDOS has a finite constant value for zero en- 
ergy (as for the monolayer) and it also shows the ZO' 
peak at 94 cm"^ (as in graphite). 

Let us consider the broadening of the LA and TA 
branches in Fig. [9] As already said, in graphene, these 
broadening do not vanish for q — >■ and, for small q, 
they are relatively constant over a wide range of q values 
{e.g. for the LA mode we are considering the region with 
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FIG. 9. (Color online) Anharmonic phonon broaden- 

ing (FWHM) at 300 K. Graphene acoustic branches (solid 
line) are compared with the corresponding graphite acousti- 
cal (dotted) and quasi-acoustical (dashed) ones. 



FIG. 10. (Color online) Anharmonic phonon broadening 
(FWHM) at 300 K. Graphene monolayer acoustic branches 
(solid line) are compared with the corresponding bilayer 
acoustical (dotted) and quasi-acoustical (dashed) ones. 



q < 0.4M and q < 0.4K in Fig. |9]). This, "plateau" is 
due to Normal scattering towards the ZA phonons and 
its characteristics stem for the fact that the ZA disper- 
sions is quadratic and that the integration (the sum in 
Eq. p| is done on a two dimensional Brillouin zone. On 
the contrary, for three dimensional graphite, the broad- 
ening of both LA and TA modes vanishes for g — >■ 0. 
It is remarkable, however, that the graphite broadening 
still presents a Normal scattering plateau, similar to the 
one of graphene, for sufficiently large q. This is par- 
ticularly evident for the LA mode for q > 0.2M and 
q > 0.2K in Fi g.[9| The LA and TA broadening in bilayer 
graphene. Fig. 10 are rather more similar to the graphite 
one than to the graphene ones, indicating that, for a 
higher number of layers, the broadening should rapidly 
converge to the bulk graphite one. Note that, graphene 
bilayer presents nonvanishing broadening of the TA and 
LA bands at F, its magnitude being about half for bilayer 
than for graphene. 

Finally, let us consider the z polarized branches. The 
ZA broadening in graphite is similar to the graphene one. 
On the other hand, the ZO' broadening of graphite is 
much larger than the ZA one, in spite of the fact that 
the ZO' and ZA branches are strictly related. Particu- 
larly striking is the sudden increase in the ZO' broaden- 
ing for certain values of q in the vicinity of F (e.gr for 
q = 0.47M along the FM direction). This peak of the 
broadening is due to the decay of a ZO' phonon, having 
a finite wavevector q into a ZA phonon near q and a LO' 
(or TO') phonon near F. This kind of decay is possible 
only when the energy difference between the ZO' and the 
ZA is equal or smaller than the energy of the LO' and 
TO' at F (the LO' and TO' are degenerate at F), see 
Fig. [TT] This condition is verified only for q sufficiently 
far from F. Thus, by increasing g, the sudden availabil- 



ity of this new decay channel produces the peak in the 
broadening. Finally, for the bilayer, the ZO' broadening 
presents a structure which is not fundamentally different 
from the one in graphite. 



C. Thermal transport 

The intrinsic anharmonic thermal conductivity kl bas 
been computed within the single mode time relaxation 
approximation using Eq. [7] . For the two dimensional 
materials (graphene monolayer and bilayer) we have used 
the convention that the volume 17 in Eq. [7] is the surface 
planar unit cell multi plie d by the inter-layer distance of 
graphite, 3.32 A. Fig.] 
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reports the thermal conductiv- 
ity and its decomposition into different branch contri- 
butions. In the temperature range considered, the con- 
ductivity is almost entirely due to the acoustic and the 
quasi-acoustic branches. Calculations of Fig. [12] arc done 
for T > 200 K. Below that temperature, converged re- 
sults can be obtained only by using a much finer q-points 
grid than those presently used. 

Let us, first, consider graphene. From Fig. [T2j the 
ZA contribution increases by decreasing the temperature 
(actually, it diverges for T — > 0), while the LA and TA 
contributions are non monotonic and reach a maximum 
near 300 K. The difference in the two behaviors can be 
understood by considering that, for small T, the phonons 
mostly occupied have small g, and that, for g — >■ 0, the 
anharmonic broadening (the inverse of the t appearing 
on the r.h.s of Eq. [?]) of the ZA mode goes to zero at 
any T. This is not the case for the broadening of the 
TA and LA branches. Fig. [3] . Now, let us compare in 
Fig. |12| graphene with graphite. The ZA contribution in 
graphene corresponds, in graphite, to the two separate 
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FIG. 11. (Color online) Phonon dispersion and anharmonic broadening of the acoustic and quasi-acoustic branches of graphene, 
graphite, and graphene bilayer. The upper panels show the lower part of the phonon dispersion as in Figs. [2] and [8] The lower 
panels show the FWHM at 300 K. The values in the A-F section for graphite are magnified by a factor 10. 
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FIG. 12. (Color online) In-plane thermal conductivity obtained within SMA (thick solid line) and its decomposition into the 
contributions due to different phonon branches. The figure shows only the contributions from a subset of the most relevant 
phonon branches. The sum of these partial contributions is the thin solid line close to the total conductivity. In order to 
compare more easily the curves, the single branch contributions of graphite and bilayer are scaled by a factor 2. Dots are 
experimental data from Ref. 1291 



contributions ZA and ZO'. These two are quite different 
already at room temperature. Below room temperature, 
the ZA increases and diverges for T — > (as for the ZA in 
graphene), while the ZO' does not. The TA contribution 
in graphene corresponds, in graphite, to the TA and TO' 
ones. Above 200 K, the TA and TO' contributions are 
almost indistinguishable, in spite of the fact that only the 
TA branch is actuaUy acoustical. Important differences 
between TA and TO' contributions appear only below 
50 K (not shown). The same considerations hold for the 



LA LO' couple. By comparing in Fig. 12 graphite with 
the bilayer, above 200 K, the overall behavior of the two 
systems is relatively similar, in spite of the different di- 
mensionality. Indeed, in the same temperature range, the 
total conductivity of two dimensional graphene mono- 



and bi-layer and that of three dimensional graphite are 
relatively very similar. Fig. |13[ with a difference between 
graphene and graphite of less than 10%. 

Fig. [12] also reports the measured in-plane thermal 
conductivity. This is done only for bulk graphite, be- 
cause of the abundance of experimental data. At present, 
experimental measurements on graphene exist only for 
small samples, where border effects are important, con- 
sequently the range of measured values is large and 
the number of samples in temperature limitecpSl_ From 
Fig. 12 the calculated graphite conductivity (which is 



obtained within the SMA) grossly underestimates the 
measured one, by about a factor two. It is unlikely that 
this disagreement is due to density functional theory. In- 
deed, DFT reproduces accurately the measured graphite 
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phonon dispersion, Fig. [T] suggesting that the most im- 
portant quantities used in Eq. [7] are correct. On the 
other hand, as aheady explained at the end of Sec. |II Al 
the thermal conductivity calculated according to Eq. [7] 
derives from the single mode relaxation time approxi- 
mation and not from an exact solution of the transport 
equation. Indeed, according to Refs . [?7| and 1501 the SMA 
cannot be used to properly describe the in-plane thermal 
conductivity in graphitic materials and the exact solution 
of the Boltzmann transport equation is required. How- 
ever, the results of Refs. [17] and [30] are obtained by us- 
ing a semi-empirical interatomic potential and a direct 
comparison with the present results is not meaningful. 
Further investigation on this point is required. 

We now consider the thermal conductivity along the 
z axis, perpendicular to the graphene planes. Fig. [M] 
shows calculations and compares them with measure- 
ments. The quasi totality of the conduction is due to 
the acoustic and quasi-acoustic phonons polarized along 
z and, as expected, the conduction is much smaller than 
the in-plane one. The transport in-plane and the one 
along z are quite different. The phonons relevant for the 
z conduction have much smaller velocity and are local- 
ized around 7 in reciprocal space. Indeed, only phonons 
with a nonzero dispersion along the z-axis can conduct 
since they have a nonzero velocity and, thus, can give a 
contribution to the sum in Eq. [7J These phonons belong 
a small region surrounding the F-A line. If we only inte- 
grate the contribution from a zz-oriented q-space cilinder 
centered around F, we estimate that, at 300 ,K the cen- 
tral 20% of the WS cell contributes about 85% of the 
transverse transport, but only aout 40% of the in-plane 
transport. 

We also remark that the the z conduction is extremely 
sensitive to the value of the cja lattice parameter. In- 
deed, a small change in cja results in a systematic in- 
crease or decrease of the frequencies of all the phonons 
relevant for the transport along z, see Fig. [T] Moreover, 
a systematic rescale of phonon frequencies by a certain 
factor A results in a rescale of the conduction by a factor 
which can be much bigger than the initial A (see Eq. Uj- 
The agreement with measurements from Fig. [M] is sat- 
isfactory; the theoretical calculations slightly overresti- 
mates the experimental value, as we expect since we are 
omitting isotopic effects. We judge it compatible with the 
assumption that the SMA correctly describes the thermal 
transport along the z axis. 



IV. CONCLUSIONS 

We have developed and implemented a generic method 
for the calculation of anharmonic three-phonon scatter- 
ing coefficients within density functional theory. The 
three phonons can have three arbitrary wavevectors 
(q, q',q"). The approach works for materials with an 
electronic gap and also for metallic or zero gap materi- 
als. The method has been implemented in the Quantum 
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FIG. 13. (Color online) In-plane thermal conductivity calcu- 
lated within the SMA for for graphene single and bilayer and 
bulk graphite. Measurements (exp.) are done on graphite 
and are from Ref. ,29-. 
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FIG. 14. (Color online) Thermal conductivity of graphite 
along the direction orthogonal to the planes and its decom- 
position into different phonon contributions. Measurements 
(exp.) are from Ref. 1291 and should be compared with the line 
labeled as "total". 



ESPRESSO package^ and generalizes a previously ex- 
isting code developed in Ref. [HI The anharmonic coef- 
ficients which can be obtained can be used in a straight- 
forward way to compute the anharmonic broadening of a 
phonon with an arbitrary wavevector q and the intrinsic 
thermal conductivity within the single-mode relaxation 
time approximation (SMA). The first applications have 
been devoted to study graphite, graphene mono- and bi- 
layer. 

We have reported a detailed analysis of the anhar- 
monic phonon broadening in graphene. Interesting, the 
broadening of the high-energy optical branches is highly 
nonuniform and presents a series of sudden steps and 
spikes of various origin. At finite temperature, the two 
linearly dispersive acoustic branches TA and LA have 
nonzero broadening for g — > 0, as noticed in Ref. jS) This 
anomalous behavior is due to Normal scattering towards 
two ZA phonons (which are quadratically dispersive), for 
small q. The activation of the Umklapp scattering for 
sufficiently large q is associated with a sudden increase 



of the broadening at nearly half the Wigner Seitz cell. 
We provide a set of expressions which ca be used to fit 
the anharmonic scattering time and broadening for the 
acoustic phonon branches, which are the most relevant 
in thermal transport. 

The broadening of graphite and bi-layer is, overall, 
very similar to the graphene one. The most remark- 
able feature is the broadening of the quasi acoustical ZO' 
branch, which is much larger and very different from the 
one of the strictly related ZA acoustic branch. On the 
other hand, the broadening of the TA and LA branches 
of graphite, displays a certain number of similarities with 
that of graphene mono and bi-layer, in spite of the dif- 
ferent dimensionality of the systems. 

Finally, we have calculated the intrinsic anharmonic 
thermal conductivity within the SMA. The in-plane con- 
ductivity in graphite, graphene mono- and bi-layer are 
very similar in spite of the differences among these sys- 
tems. The calculated SMA conductivity heavily under- 
estimates the measured one (for graphite) by almost a 
factor two, in the temperature range from 200 to 2000 
K. This finding is compatible with the conclusions of 
Ref. [57| and 1301 which state that the SMA cannot be 
used to properly describe the in-plane thermal conduc- 
tivity in graphitic materials. On the other hand, the 
calculated SMA conductivity for graphite along the di- 
rection perpendicular to the plane is in good agreement 
with measurements. 
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Appendix A: Third order calculation 

This section describes the method used to calculate the 
third order anharmonic scattering coefficients. In order 
to fix the notation, Sec. |A1[ and A2 resume DFT and 
linear response to DFT (alias DFPT). Third order cal- 
culations are described in Sec. |A3| Sec. |A4[ and |A5 
gives the explicit expressions for certain terms. Sec 



X6 



describes the implementation of the nonlinear core cor- 
rections. 



1. Kohn-Sham equations 

Within DFT the total energy a system can be deter- 
mined from the ground state electronic charge density n. 
In turn, n can be obtained by solving self-consistently 
the Kohn Sham (KS) equations'^-'^ , which, in a periodic 



crystal are: 



T/'^==(r)=w'™(r) + 



5Ei\n] 



5n{v) 
i(r) = ^0\,|(V'k,z|r)P; ( n{v)dv 



N' 



el 
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(Al) 
(A2) 

(A3) 



In Eq. Al r"*'" is the single-particle kinetic energy oper- 
ator, V^^ is the KS potential, iV'k.i) are the Bloch eigen- 
states with wavevector k, band index i, and energy ek.i- 



(r + RlV-k.^ 



ik-R, 



r|'0k.i)j being r the position and 



R a lattice vector. u'°" is the external potential due to 
the ions, Ei[n] is the interaction functional (Hartree en- 
ergy plus exchange-correlation contribution). The sum 
in Eq. |A3| is done on a sufficiently fine grid of k-points. 
0k, i is the occupation of an electronic state: 0k, i — 1 for 
valence band electrons and 0k,i = for conduction ones. 
Here and in the following J dv is the integral over all the 
space. iV' is the total number of valence electrons (we 
consider e^ — 1). The total energy of the system is: 



£;tot ^ ^1, 



k,i 



<^-k,iO\^A+E\[n]- 



5Ei[n] 
(5n(r) 



n{Y)dr, (A4) 



where i^'°" is the ionic contribution. 

In this form, the KS equations are suitable for semi- 
conductor or insulators, where the electronic gap is differ- 
ent from zero, if the electronic gap vanishes (metal and 
semi-metal case) it is customary22l ^ introduce a smear- 
ing function 6„{x), which is characterized by a smearing 
width a^ and which becomes a step function in the limit 



a — » 0. The KS equation are still written as in Eqs. ~Kl 
A3 but now 0k, j = Qa{^F — ek,i), where the Fermi energy 
ep has to be determined self consistently from 



E' 



'k,j 



iV° 



(A5) 



Furthermore, in the metallic case , a proper definition of 
the energy £^'°* requireP^l Eq A4 to include the term 



E 



where 5a{x) — d9„{x)/{dx). 



xSa-{x) dx, 



(A6) 



2. Linear response 



The derivative of the electronic charge distribution 
with respect to the q periodic displacement, Uq as de- 
fined in Eq. nl (for simplicity from now on we will drop 
the indexes a and s}, can be obtained from first order 
perturbation theoryI21, Por the metallic case, linear re- 
sponse can be implemented following Ref. 1331 we will 
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follow the equivalent, but slightly different approach of 
Ref-lMi 

Let us consider a uniform grid of electronic k points 
and a phonon wavevector q which belongs to this grid. 
First, one has to solve the KS equations and obtain the 
ground state change density n and the corresponding KS 
wavefunctions |V'k,i)- Then, one has to define a "cutoff" 
energy E which separates the electronic states which are 
completely empty from those which are occupied or par- 
tially occupied. In the semiconductor/insulator case, E 
can be any energy within the gap; in the metallic case 
any E > ep + 3a is a, reasonable choice. We define Pc 
as the projector on the manifold spanned by the empty 
states and P^ = 1 — Pc as the projector onto the occupied 
and partially occupied states. 

The derivative of the charge dn/duq and the derivative 
of the KS wavefunctions projected onto the conduction 
bands j^^ -) = Pcldipk^i / du^) can be obtained by solving 
self-consistently the equations: 



[T' 



kin I T/KS 



aPv - ek,i]|</'kj) = ^^c-^ IV'ki 



dur, 



(A7) 



dUq 



9i;'°"(r) 



5'^Ei[n] 
(5n(r)(5n(r' 



9n(r') 
duq 



dr' (A8) 



dn{r) 
dua 



^('^\{ X!'^k»eFl'/'ki)(V'kJ 



7ki ^ "'k+qj" 



+ E T. T''''. IV^k+q,J-)(V'k+qJV-q|Vk.)(V'k.| 



«J 



Ckj ~ ^k+qj 



E^~k.[iO(V'k,i + iv^k.)(0k;'i 



(A9) 



a is a const ant chosen in such a way that the linear 
system of Eq. A7 is not singular. ^^ indicates that the 
sum is to be performed only on the partially occupied 
states. ^k,i = ^<T(eF - ek^i) 

The first two terms in the right hand side of Eq. X9 
are different from zero only in the metallic case and are 
written using the notation Cp = dep/duq and V^ = 



dVKs/du^ 



,q _ 



for Qt^O; it has to be determined 



self-consistently from 



Ek.,:(^k,»|V^'^|V'k, 



(AlO) 



In the second line of Eq. |A9| we have used a compact 
notation which reads: when the denominator is equal to 
zero one has to substitute the ratio with its limit as the 
denominator approaches zero. The same substitution can 
be used when the denominator is very small in order to 
gain numerical stability. Thus, when ekj ~ ^k+qj" one 
can substitute (6'kj - 6'k+q,j)/(eki - Ek+qj) with -Si^j. 
We, finally, note that the present approach is different 
from the the one described in Ref. [551 and that the \(j>) 
wavefunctions presently defined are different from the 1 0) 
of Ref. |33l 

3. Third order 



Let us consider three phonon displacements Uq, Uq', 
Uq" such that the sum of their wavevectors is q -f q' -I- 
q" == 0. By solving the linear response equations one 
can obtain dn/duq and the {10^ ^)} corresponding to the 
three phonons. The third derivative of the energy can 
then be obtained from 



d^£' 



dUqdUqidUqn 6 



^q,q .q ^ ^q .q.q _|_ £;q .q .q _|_ ^^q-q -q ^ £;q .q -q _|_ £;q ^q.q 



(All) 



^q.' 



^q,q 



^3^ion 
dUqdUq'dUqii 



dn{r) a2t>'°"(r) 

O / —T\ -7. 7{ "I" 

aUq OUq'OUq" 



"W75 S A '^^ 

OUqOUqiOUqii 

5^Ei[n] 



dn{r) dn{r') dn{r"] 



(A12) 



Sn{r)Sn{r')Sn{r") duq duqi duq 



dvdv'dv" 



For the semiconductor/insulator case, we can follow Ref. [13] and write: 



2'qq q 



k 



Y.~^U4>i^\v'^\4l) - E^~k.(c^q,i0tq,,)(v^k-q',iW|^k, 



*j 



(A13) 



For the metallic case, we follow Ref. [HI 
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k I i 



'M{c^^?\V''\'t^li. 



6E 

V 

2E 



^ [4+q,.(C+q,.l^''>k-q'.,) -e~k-q',,(^k+q..|^'>|0lq,,,)] (^k-q' J I ^-l" |^k+q,.) 



(^k,^|^''|V'k^qj)(^k-qj|l^''>k+q",/)(V'k+q",i|l^''"|V'k.,:) 

^k,i(ek-qj — fk+q",;) + ^k-qj(ek+q",/ — ^ks) + ^k+q",i(ek,i — fk-qj) 

(ek,i — ek-q,j)(fk-qj- — ek+q",/)(ek+q",; — Ck.i) 



3e^ 



-^ <Jk,i ^ Ok+q" J 



E 



ek,i — Ck^ 



^(^k.d^'^1V'kw^J■)(V'k+q'^,|W'|^k.) + 2^4,»(^k,.|^'''|<-: 



344 ES^^mI^'^'^m) -e^e^'ef K] ^, 



(1) 

k,i 



(A14) 



where 4^) = 9(5^(2:)/(9a;)|^^,p_,^^ 



Eq. 



A14 



is written with the same compact notation of 



A9 



when one of the 



denominators in Eq. A14 vanishes the corresponding term is replaced with its limit. In particular, when Ck 



£k-q',j hi the second line of Eq. A14 the argument of the sum can be written as 



q,' 



5'k+q,*(</'k+q.»l<^k-q'j) 



~Sk+^.A^k+^.^\V''K-^',,) ('Ak-q',,l^'''>k+q,,) 



The limits of the factor in the fourth line of Eq. |A14| are 

for ek,i ^ £k-q,i 7^ £k+q",; : 
9k, i — ^k+q",; J 



ek,i ~ £k+q",i 



1 

^k+q",; ^ Ck,! 



fk,; 



for Ck-qj ^ Ck+q",/ ¥" £k,i : 

<5k-qj 

for Ck+q",; ^ ek,i 7^ <:k-qj : 

^k+q",i 



t'k-q,j — t'k,! 
^k— q.j ^k,7 



'k+q",i ^ '^k-q,j 



Ek-c 



fk^ 



,i — Ek-c 



1 



£k-q,j ^ Ck+q",; 



for ek.i 



ek-q,i "^ ek+q"j 



1?(1) 

2^,.- 



(A15) 



Finally, in the fifth line of Eq. A14 when eki ^ ^k+q".j 
one can substitute ((5ki — Sk+q".j)/i^'ki — ek+q".j) with 



-m- 



Once the the derivative in Eq. |A11| has been deter- 
mined, one can obtain the phonon scattering coefficients 
by combining Eq. |4]and Eq.[5| Provided that G is vector 



I 

of the reciprocal lattice, we also remark that 

Q3ctot nSctot 

dUadUa'dUai'+G dUqdUqi dUq" 



(A16) 



Hence we have not lost of generality by imposing q + q' + 
q" = at the beginning of the present section. 

Given a computer code which implements linear re- 
sponse to DFT (DEFT), all the bra-ket products de- 
scribed in this section can be obtained straightforward. 
On the other hand, the computation of the second to the 



fourth terms in the r.h.s of Eq. A12 have to be imple- 
mented from scratch. The implementation of the fourth 
term is trivial within the local density approximation. 
The expressions required for the other terms are given 
below. 



Ionic contribution 



The second term in the r.h.s of Eq. |A12| is the third 
derivative of the ion-ion contribution to the total energy. 
It is computed, as customary, using the Ewald sum tech- 



niqui 
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Ss'.s"Zs'ZsFa^l3,j{c[,ts' — ts) + 6s" ,sZ s" Z s' Fa^p^^{ci , t^" — t^') 
+ 5s^a>ZsZs"Fa^f}^y{q[' ,ts — t^") — 5s^s',s"Zs 2_^ ^s^a,^,7 (0, t^ — t^). 



(A17) 



In Eq. A17 we have written explicitly the dependence 
on the atomic {s,s' ,s") and Cartesian {a,j3,^) inde xes of 
the phonon patterns Wq (defined in Eq. II]). In Eq. A17 



Zg is the ionic charge and t^ is the position of atom s. 
The sum is performed over all the atoms of the unit cell. 
The function F is 



^a,/3,7(q,t) = - 



G 



2-(G+q)V(V 



i(G+q)-t 



Xz'''(G + q)„(G + q)/3(G + q), 



R 



,iqR 



d3/(x) 



dxadxpdx-y 



(A18) 



x=t-R 



Here, fl is the unit-cell volume, the sums are performed 
on the ensemble of the reciprocal lattice vectors G and 
of the real space lattice vectors R. r] is the cutoff for the 
real space summation within the Ewald methocP^ and 
/(x) = erfc(?7|x|)/|x|, being erfc is the error function. 
The derivative of / is 



rfV(x) 
dxadxpdx^ 



{6apXj + Sfj^Xa + S^aXj3)fi{\yi\) 



+ XaXi3Xjf2{\:X.\) 



(A19) 



with 



f2{x) 



3erfc(?7a;) + a(?7a;)(3 + 2ri x ) 
x^ 
15erfc(77x) + 0(770;) (15 + lOrj'^x^ + irj'^x^) 

x^ ' 

2^ .2 



5. Derivatives of the external potential 



Here we give the expressions to calculate the third and 



fourth terms in the r.h.s of Eq. A12 Both terms are 
convenient to evaluate in the reciprocal space. w'°'^ in 
Eq. |A12| corresponds to the electrostatic potential in- 
duced by the atomic ions. Within the present approach 
an atom s acts through a pseudopotential, which, within 
the Kleinman-Bylander scheme '^, is nonlocal and can 
be written 

(A20) 



where v°'^ is the local component of the potential; D^ ^ ~ 
Df, are coefficients and P^^ are ion-centered projectors. 
The total ionic potential is a superposition of ionic po- 
tentials: 



w'°"(r,r') 



R,s 



R,r' 



R) 



(A21) 



where the sum is performed on all the lattice vec- 
tors R and on the atoms s in the unit cell. When 
w'°'^ is local, its trace with the charge density is 
J u'°"(r)n(r)(ir. When w'°" is nonlocal, the same quan- 
tity is ^^k,Jdr/drV*^(r)v'°"(r,r')V^;:,(r'). With 

this notation 



/ n(r)t;'°"(r)dr = i ^ n{^GWr{G)e 



-iG-ts 



G,s 



A. 






V';;,,(G)P^,.(k + G)e 



— JG-t., 



(A22) 



(A23) 



Here, n.(k), w^°^(k), P^_s(k) are the Fourier transform 
of ?T.(r), wl,°'^(r), P^ s(r). The Fourier transform of /(r) 
is defined as /(k) = J e~'^^''^f{r)dr, where the integral 
is done all over the space. Given a Bloch wavefunction 
V'k,i(r), we define Vk,»(G) = / e-*(''+<=^>''V'k,i(r)dr, that 
is Vk.,(r) = l/{Nn) J: e^(k+G).r^^^^(G). 

G 

To simplify the notation we notice that the derivative 



of the charge, Eq. X9 can be rewritten as 



|^(r)=5:0M{0S,(r)Vk,(r)-fVk.(r)fe,?(rr}, 



and we define 



^UG) 



,-l(k+q+G)-r7q 



^U^)dr. 



By using the above definitions, the third and the fourth 
terms in the r.h.s of Eq. |A12| are 



/ 



n(r) 



gabion ^j.) 



^Uq^s,a^Ucl'^s\|3^Uc^",s",'y 



-dr = 



^s,s' ,s" 

n 



Y,ni-G)i-ifG^GpG^vl°%G) 



„-.G.t. , ^-^'-"J2d;J, 



G 



X iA 



(",/3,7) 

fJ,,S 



A, 



A 



(",/3) 
k,i 



A 



(7) 



A 



(a) 

k,i 

/J.,S 



A 



(/9,7) 



N 






A 



(a, 0,7) 
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(A24) 



dn{r) d^v'°'\r) 

dUq^s,a t^Wq'.s'./jSMq",^" 



-dr 



E ?^^-^(-^)'(q + G),(q + G).^^i'?^(q + G)e-(^+«)-*=' +2^fY: 



^ ^ ^ *-'^— q,s,Q 




B, 



q.k,! 



A 



m 



B 



if) 

q,k,! 



A 



(7) 



B 



(/3) 
q,k,! 



A, 



^l,s' 



B 



(/3,7) 

q,k,i 



/i,I/ 



(A25) 



In Eq. A24 the overline indicates the sum on the per- 



mutations of the Cartesian indexes, e.g. A°'l^[Ai]* — 

^a/3[^7]*+^/37[yi"]*+^7"[^/3]*. Furthermore, ^^ is 

the Fourier transform of » ■ The A and B coefficients 



of Eqs. A24 and A25 are defined as 



A 



(a) ^ -I 

k.i 






n 



E^L(G)e-«-*= 



G 



x(k + G)„(k + G)0P^,,(k + G), 

,(a,/3,7) _ (-i)^ 



■^k,i 

IJ,,S 



n 



E^L(G)e-«-*' 



G 



X (k + G)„(k + G)/3(k + G)^P^,,(k + G) 



G 



i?l'!h = ^Efe(Gr(k + q+G)„ 
G 

xF^^,(k + q + G)e-*('i+°)-*% 



B 



(c/J) _ J-i)' 



l.k.i 



r2 



Efe(G)]*(k + q + G)„(k + q+G) 



G 



xP^,,(k + q + G)e-*(i+°)-*=. 

Where the notation (k + G)q means: The cartesian com- 
ponent a of vector k + G. 



6. Nonlinear core correction 

Within the pseudopotential approach, the electronic 
charge density is divided into core and valence contri- 



r 



butions. A common way to include the core effects in 
the Kohn-Sham equations is the nonlinear core correc- 
tion scheme of Ref. 1361 In practice, one determines the 
charge density of the core electrons nc(r) only once, be- 
fore the KS self-consistent cycle, as nc(r) only depends on 
the pseudopotentials. The KS equations are still solv ed 

re- 



only for the valence bands and n(r) defined in Eq. A3 
mains the valence charge. On the contrary, in Eq. |A2|one 



substitutes Ei[n\ with £'i[nt] where nt{r) — n(r) + nc(r) 
is the total charge density. 

By using this approach the third order equations need 
to be modified. First, in the last term in the right hand 
side of Eq. A12 one has to replace the three dn/du^ with 



the corresponding dn^/duq. Second, one has to include 
in the r.h.s. of Eq. |A12| the two additional terms 



-dr 



Sn{r) dUqdUqidUq'i 

g r S^E,[n^] d'n,{r) d^n.jr') ^^^^, 
J Sn{r)Sn{r') duqduq> duqi> 



(A26) 



Appendix B: Notes on the implementation 

This section discusses some practical issues concerning 
the actual implementation of the method. 



1. Implementation 

Let us consider three wavevectors in the general case in 
which q ^ q' ^ q", with q+q'+q" = 0. The calculations 
of the third order derivatives is done in three consecutive 
steps. 

1. Self consistent calculation to obtain the ground 
state charge n(r) and the Kohn-Sham potential 
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2. Self consistent linear response calculation to deter- 
mine dn{r)/duq and dV^^{r)/duq. This is done 
three times for the three cases q = q', q", and q". 

3. Calculation of the third order derivatives. 



In order to do Step 3, seven distinct sets of ground 
state wavefunctions lip) are needed (see eq. A13 or 

-q ^ 



A14I 



IV'k), IV'k+q), l^k-q), with q=q, q', and q", where k 
runs on the ensemble of wavevectors which are used in 
the electronic integration. The jf/'k) are already calcu- 
lated in Step 1. They can be saved and read from disk or 
they can be calculated a second time with a relatively in- 
expensive non self-consistent solution of the KS equation 
(based on the knowledge of V^^ and 7i(r) determined in 
Step 1). When k runs on a regular grid and q belongs 
to the grid, the different sets might coincide. However, 
given the way the code is implemented U-L (the wavefunc- 
tions are read and written from sequential files; different 
processors might work with different k points in parallel) 
it is bests to keep distinct the seven sets. 

Twelve distiiict sets_ of wavefunctions derivatives are 
also needed: |0^), Icj)'"^), |0^_q), and It^-k+q). with q=q, 
q', and q". These wavefunctions can be calculated with a 
non self-consistent solution of the linearized KS equation 
which is based on the knowledge of the dV^^{r)/duq cal- 
culated in Step 2. Even in this case, for practical reasons, 
it is better to keep distinct the twelve sets even when 
some of them coincide. These derivatives are computed 
resolving non-self-consistently eq. |A1[ in the present 
work, this was the most CPU intensive step; on the other 
hand, for bigger systems the Input/Output of dn/duq 
can become a bottleneck. 

In order to evaluate Eq. |A13[ one needs to calcu- 



late matrix elements of the kind 



lyq 



and 



(V'k-q'l^'^ IV'k+q)- In both cases this calculation has 
to be done six times for all the possible permutations 
of the three wavevectors q, q', q". For the metallic 
case, Eq. A14 further terms are needed. The terms 
(^ic+q l^'^'IV'k+qt) need to be computed six times for 
c[a and qb equal to q,q', q", with qa ^ q?,. The terms 
(V'k+q|l^''|V'k) and (^k-q|^~''|V''k) are computed already 
in Step 2 and can be saved for later use in Step 3. The 
terms {ipiti-qlV'^ jf/'k+q") need to be computed for the six 
permutations of q, q', and q". 

In some special cases the number of operations can 
be greatly reduced. The most obvious one is when 
q = q' = q" = 0: in this case only one set of wave- 
functions derivatives is needed: |</'k)- Moreover, when 
q = 0, q' = — q" = p, with p ^ 0, only three sets 
of wavefunctions derivatives are needed: |0^), \(j)^}, and 
^k-q)- This is because 9n(r)/9M_q = [dn{r)/duq\* and. 



:^(r) = [0^(r) 



I'/'k-q)- -^"iC"" Ljcv.a,uoctyin_±;/ty«_q 

according to time reversal symmetry. 
Another special case is when q = 2p, q' = q" = — p; in 
this case 8 derivatives are needed and index permutation 
can be used to avoid computing some terms. 



2. Symmetry and q-points 

A Bravais lattice can be invariant under a certain num- 
ber, up to 48, of rotations. A crystal will respect a sub- 
set, possibly all, of these symmetries defining the group 
of crystal symmetries, Q. These symmetries can be ex- 
ploited to reduce the computational cost of the calcula- 
tion and must be enforced to avoid unphysical breaking 
of symmetry caused by computational noise. Here, we 
briefly revise how this is done in the pwscf and phonon 
codes of Quantum ESPRESSO, and we describe how 
the approach has been extended to the third order. 

In the ground-state energy calculation, the Kohn-Sham 
equations (Eqs. |Al[ A2 A3), are not solved for all the 



k-points of the chosen grid. Instead they are solved only 
for a subset of k-points which are inequivalent under the 
rotations of G (these k points belong to the so called irre- 
ducible wedge of the Brillouin zone) . The charge density 
n(r) is then obtained, inexpensively, by summing on this 
k-point subset and, then, by imposing the symmetries of 
the crystal. The symmetries are imposed by rotating the 
initial charge density according to all the symmetries and, 
then, by making the average. In general, even if the crys- 
tal has no symmetries, the sum can be performed on half 
of the k-points of the initial grid since ■0_k(r) — ''/'k(i") 
(time reversal symmetry). 

In the phonon calculation, the situation depends on 
the phonon wave- vector q. Let us consider the the sym- 
metries of the crystal which are associated with a rota- 
tion that leaves q unchanged (or that transform q into an 
equivalent q + G). This subgroup of Q is called the small 
gro u p of q, Gq- The linear response equations (Eqs. A7 
AS A9) are not solved on the grid. Instead, they are 



solved only for a subset of k-points which are inequiva- 
lent under the rotations of Gq- dn{r)/duq is obtained by 
summing on this k-point subset and, then, by imposing 
the symmetries of Gq (in analogy with n(r) in the pre- 
vious paragraph). Because of time reversal symmetry, 
in order to obtain dn{r)/duq from Eq. A9 one needs to 
compute explicitly only the |(/>^) but not the \4>^'^)- How- 
ever, if there are no other symmetries, the sum has to be 
performed on all the k-points of the grid and not on half 
as in the total energy calculation. Once that i9n(r)/9uq 
has been calculated, one can determine the dynamical 
matrix at q. 

Another way in which the symmetries can be exploited 
is the following. The ensemble of the vectors obtained 
by rotating q with all the rotations of G is called the 
star of q. The number of inequivalent vectors in the 
star is |5|/|5q|, where \G\ is the order of G- Once that 
dn{r)/duq has been determined, all the dn{r)/duq for 
the vectors q in the star can be obtained inexpensively by 
rotation. Moreover, because of time reversal symmetry 
dn{v) / du-q = [9n(r)/(3Mq]*. 

A third-order calculation is conceptually not differ- 
ent. In this case we are dealing with a triplet of points 
(q, q', q") and one has to consider the small group of the 
rotations of G that leave the three q-points unchanged. 
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laicr 



The electronic k-point summation in Eq. A13 



A14 (and also those in Eqs. A24 A25 1 is done on a 



subset of point of the initial grid which are inequivalent 
under the rotations of ^q,q',q". The third order matrix 
obtained from this partial summation is not useful as it 
is. The actual third order matrix is obtained by imposing 
on it the symmetries of ^q.q'.q"- 

In analogy to the dynamical matrix case, once the 
third-order matrix of a given (q, q',q") triplet has been 
determined, on can apply the rotations TZ of the crystal 
symmetry group Q and obtain inexpensively the matri- 
ces corresponding to the rotated triplet {TZc[,TZq' ,TZq"). 
Moreover, the matrix of the triplet (— q, — q',— q") can 
be obtained by conjugation. Another useful symmetry 
is the trivial one associated to the index permutation: 
by construction, the matrix associated with the triplet 
(q, q',q") is equal to the six matrices obtained by per- 
muting the three vectors q,q',q" and the corresponding 
indices. 



The triplet obtained in this way are deleted from the list 
of the "triplets to be done" and will not be computed in 
the following step of the loop. This procedure allows for 
a spectacular reduction in the number of triplets; in the 
graphene case 4096 possible triplets are reduced to just 
88 independent triplets. 



Appendix C: Fourier interpolation 

Actual DFPT calculations are done on a relatively 
coarse grid of q wavevectors. Dynamical matrices and 
third order coefficients are then obtained for a finer grid 
with a Fourier interpolation technique. In this section, 
first, we revise the Fourier interpolation technique as it 
is implemented in the standard Quantum ESPRESSO 
package. Then, we describe how the method has been 
generalized to third order force constants. 



3. Discretisation on a grid 

In some cases, it is useful to calculate the third-order 
matrices on a regular grid of q wavevectors in the Bril- 
louin zone, the reason for this are clarified in Sec. [C] 
When we say that calculations are done on a given grid, 
it means that we have calculated the third order coeffi- 
cients corresponding to every triplet (q, q', q") of points, 
such that each vector belongs to the grid and that the 
condition q + q' + q" = G is satisfied. 

In the actual implementation of the procedure, we first 
determine dn(r)/duq and the dynamical matrices on the 
grid of q wavevector. Exploiting the symmetries can al- 
low an important reduction of the computational cost. 
Indeed, the actual ab-initio calculation is done only for 
those q points which do not belong to the same star (de- 
fined in the previous section), or which are not equivalent 
to the opposite of a point already calculated. dn{r)/duq 
and the dynamical matrices for the other points are then 
obtained by rotation or conjugation. Once this is done 
we perform the third order calculations. 

In practice, we can perform a double loop on the grid 
and chose q and q' so that they belong to the grid. The 
third vector is chosen as q" = — q — q'. q" may not 
belong to the grid, q" is however still connected to a 
point in grid by a reciprocal lattice vector G and, thanks 
to Eq. |A16[ we are not losing generality. Moreover, 



dn{r) 

dUq" 



dn(r) 



d 



q"+G 



for every reciprocal space vector G, and we can use the 
i9n(r)/9uq calculated on the original grid. Once that a 
triplet has been computed we can determine for free all 
the triplets which are equivalent by rotation, permuta- 
tion of the indices, or by conjugation. Actually, some of 
these operations can be redundant, e.g. a certain rotation 
could be equivalent to a permutation, or to a conjugation. 



1. Second order 

Let us consider a lattice with basis ai, a2, a3. The 
dynamical matrices D2 ( ^ ^' ) (we use the definition of 
Eq. [2] and we drop for simplicity the Cartesian indexes) 
are first computed ab initio on a uniform grid, centered in 
the origin, of Ni x N2 x N3 q points of the Brillouin zone. 
We want to determine the real space force constants 



F. 



. s s' ^ 



d^S* 



dvo,sdvn., 



(CI) 



by Fourier interpolation of the D2 from the NiX N2X N^ 
grid. D2 for a generic q point can then be obtained by 
back Fourier interpolation. 
Let us consider 



?) = 



-E 



D2{,:>)e 



p-iq R 



(C2) 



where R is a lattice vector, A't — N1N2N3, and the sum 
is performed on the points of the grid. The F2 defined in 
this way cannot be used as they are since they are un- 

physical long-ranged constants. Indeed, F2 ( ^+^ j == 

F2 ( ^) for any R vector of the super lattice (SL) gen- 
erated by the lattice vectors Ai = A^iai, A2 = N2aL2, 
A3 = iVgas. 

The short-ranged F2 can be obtained from the F2 
in the following way. We define two vectors as "SL 
equivalent" when their difference is a vector of the su- 
per lattice defined by Ai, A2, A3 and we call Wsl 
as the Wigner-Seitz of this super lattice. Eq. |C1| is 
the force constant between the atoms whose distance is 
d = R -f tg' — ts, where t^ is the position of the atom 
s in the unit cell. We distinguish three cases: i) when 
d e Wsh, ^2 (^ ^) = -^2 (3 /); ii) when d lies on the 
border of W^sl/A ( , ?) = '^2 (, ?) /^eq, where N,^ is 
the number of points on the border of M^sl which are "SL 



equivalent" to d; iii) when d ^ Wsh, -^2 (^ ^) =0 
interpolated dynamical matrix for a generic q is then 



The 



^2(..^0=E^2(.?)^ 



iq-R 



(C3) 



R 



where the sum is done on all the lattice vectors R. This is 
the Fourier interpolation technique as it is implemented 
in the standard Quantum ESPRESSO package. 



Third order 
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Appendix D: Computational details 

1. Electronic integration 

The electronic integration of the density functional 
theory calculations is done using a first-order Methfessel- 
Paxton smearing22l of 0.02 Ry which converges for a grid 
of 32 X 32 X 1 electronic k-points in simple and bilayer 
graphene and for a grid of 32 x 32 x 8 k-points for graphite. 



For the third order coefficients the situation is analo- 
gous, although less straightforward. We use the notation: 

Q3£-tot 



D- 



( q q' q" ) 

\ s s' s" ) 



1 



N dUq^sdUq>^s'dUq"^ 



Tp / R R' R" "i _ 



d^£' 



(C4) 



dv-R,.dv-R',,'dv 



•JR" 



The matrices D^ are first computed ab initio on a 
uniform grid of q points centered in the origin, mean- 
ing that both vectors q' and q" run on the grid, while 
q = — q' — q". We then define 



FofOR'R'/N^ 



n2 2^ 

q',q" 



D- 



/ ' " 
I q q q 

\ s s' s" 



,-i(q'-R'+q"-R") 



(C5) 
where the sums are performed on the grid points. 

The force constants Fr^ of Eq. |C4| correspond to three 
atoms at the positions t^, t^/ + R', and ts" -I- R" (be- 
cause of the crystal translational symmetry we can con- 
sider without loss of generality R=0). The three atoms 
form a triangle with perimeter P = |di| + |d2| -I- jdaj, 
where di, d2, ds are the distances among the three 
atoms (the three sides of the triangle). To determine 
F3 we consider three cases, i) When all the three dis- 
tances di,d2,d3 are inside (and not on the border) of 
WsL, we put F3 ( ^; 5 ) = F3 ( ^; 5 ) . ii) If two of 
the distances {e.g. di,d2) are inside Wsl and the third 
one is outside or on the border of Wsl we calculate the 
perimeter P of the triangle. If P is the shortest perimeter 
among the perimeters of all the triangles formed by the 
triplet of atoms "SL equivalent" to the original three, we 
put F3 ( ^; 5 ) - F3 ( ° ^/ 5 ) /^eq, where iVeq is the 
number of triangles with perimeter equal to P. iii) In all 
the other cases F3 ( ° ^,' ^" ) = 0. 

The interpolated D^ matrices for a generic triplet of 
wavevectors (such that q + q' + q" = G) is obtained 
from 



D. 



E 

R',R" 



F. 



R 



L' R" ^ 

' s" J 



„i(q'-R'+q"-R") 



(C6) 
where the sums run on all the lattice vectors. The cri- 
terion presently described to determine the F^s, stems 
from the assumption of a long range exponential decay 
of the force constants. 



2. Linear response calculations 

DFT dynamical matrices are corrected using a proce- 
dure based on DFT-I-GW renormalization of the electron- 
phonon interaction as in Ref. |211 Indeed, DFT repro- 
duces very well the measured dispersions of graphite 
for all the phonon branches but for the TO one, in 
the vicinity of the high symmetry point K. This fail- 
ure of DFT, which is very specific to the graphene and 
graphite systems, has been analyzed in Ref. [Ml To im- 
prove the accuracy of the TO phonon branch, we have 
applied an electron-phonon self-interaction as described 
in Ref. ^24i The detailed procedure is described in Sec. 
IIB of Ref. [37] (third paragraph) . We have used the pa- 
rameter r = 1.65, which is appropriate to the present 
LDA calculations and which can be derived from Ta- 
ble I of Ref. 1211 Using this approach we determined 
the dynamical matrices of graphene on a super-sampled 
48 X 48 X 1 q-point grid. The matrices for the graphite 
and for the bilayer are then obtained by using as in-plane 
force constants those of graphene and as out-of-plane 
force constants those coming from independent DFPT 
calculations on the two systems. 

The third order coefficients are obtained in the stan- 
dard way, that is without including this self-interaction 
correction of Ref. [531 For graphene, the third-order co- 
efficients are calculated on a 8 x 8 x 1 q-point grid (see 
Sec. B3), which consists in 88 irreducible triplets. For 



bulk graphite, third-order coefficients are calculated on a 
8x8x2 grid, which consists in 297 irreducible q triplets. 
When computing the linewidth, we have tested conver- 
gence starting from the 8x8x2 grid coefficients, finding 
that the use of the those from the 4x4x2 subset grid (33 
inequivalent triplets) does not worsen accuracy. For bi- 
layer graphene, the third-order coefficients are calculated 
on a 4 X 4 X 1 q-point grid (12 irreducible triplets). 



3. Broadening calculations 

Eqs. [6] and [7] are evaluated by performing the sum 
over a discrete uniform grid of q points randomly shifted 
from the origin. The d{x) distribution is substituted with 
the Gaussian function S{x) = e~^^^^^ /{xV^)i where x 
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is an artificial smearing, independent from q. Tfie re- 



sults of Sects. Ill A and IIIB are obtained by using: for 
graphene, a 1800 x 1800 x 1 grid and x — ^ cm^^; for 
graphite, 600 x 600 x 15 grid and x — ^ cm~^; for the 
bilayer, a 120 x 12 00 x 1 grid and x = 2 cm"!. The 
results of Sec. |IIIC| are obtained using: for graphene, a 



128 X 128 X 1 grid and x = 10 cm~^; for graphite and 
bilayer, a 64 x 64 x 4 grid and x = 10 cm"-"^. For each 
system, the same grid is used to determine the broaden- 
ing from Eq. [6] and the thermal conductivity from Eq. [7] 
The convergence, has been tested using smaller smearing 
values and finer grids at selected temperatures. 
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